library(haven)
library(rdrobust)
library(rddensity)
library(ggplot2)
library(ggeffects)
library(dplyr)
library(scales)
library(stargazer)
library(fastDummies)
library(sjPlot)
library(fixest)
library(sjmisc)
#Load necessary packages

trial<- read_dta("trial_replication.dta") #Load the full data
civilian<-trial[trial$category_e==0,] #The data including only civilian defendants

### Table A.1. Descriptive statistics------------
t_a1<-data.frame(civilian[,c("timediff","d2h15","age","male","islander",
                             "d2_per","subversion","leakmilitary","weapon","presreyearr_l")])
stargazer(t_a1)


### Table A.2. Testing for discontinuities in covariates and types of offenses------------
m1<-rdrobust(civilian$age, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #age
m2<-rdrobust(civilian$male, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #male
m3<-rdrobust(civilian$islander, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #islander
m4<-rdrobust(civilian$d2_per, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #N of codefendants
m5<-rdrobust(civilian$weapon, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #weapon
m6<-rdrobust(civilian$subversion, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #committing subversion
m7<-rdrobust(civilian$leakmilitary, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #leaking military intelligence
m8<-rdrobust(civilian$presreyearr_l, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #president review t-1
m9<-rdrobust(civilian$d2_maj, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #major offenses
m10<-rdrobust(civilian$d2_mid, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #medium offenses
m11<-rdrobust(civilian$d2_819, civilian$timediff, kernel = "tri",bwselect="cerrd", p = 1,masspoints = "adjust") #minor offenses

r<-NULL
for(i in 1:11){
  r<-rbind(r,cbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],4)," ,",round(get(paste0("m", i))$ci[3,2],4),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
colnames(r)<-c("coefficient","robust CI", "bw","n")
r
## density test ------------------------------------------------------------
summary(rddensity(civilian$timediff))


### Table A.3. Regression discontinuity of the 1956 reform on sentence severity for civilian dissidents------------
d2h15<-civilian$d2h15
timediff<-civilian$timediff
Z<-civilian[,c("age","male","islander","d2_per","weapon","subversion",
               "leakmilitary","presreyearr_l")]
#DV, running variable, and covariates of the RD models

m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw 
m2<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m3<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw 
m4<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust",bwcheck=280) #quadratic fit, uniform kernel, CE bw
m5<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw 
m6<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m7<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw 
m8<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw 

r<-NULL
for(i in 1:8){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],3),", ",round(get(paste0("m", i))$ci[2,2],3),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],3),", ",round(get(paste0("m", i))$ci[3,2],3),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.4. Regression discontinuity of the 1956 reform on sentence severity for civilian dissidents across different types of offenses------------
d2h15<-civilian$d2h15[civilian$d2_maj==1]
timediff<-civilian$timediff[civilian$d2_maj==1]
Z<-civilian[civilian$d2_maj==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for major offenses
m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m2<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2h15<-civilian$d2h15[civilian$weapon==1|civilian$subversion==1|civilian$leakmilitary==1]
timediff<-civilian$timediff[civilian$weapon==1|civilian$subversion==1|civilian$leakmilitary==1]
Z<-civilian[civilian$weapon==1|civilian$subversion==1|civilian$leakmilitary==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for major offenses (based on crime description)
m3<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m4<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2h15<-civilian$d2h15[civilian$d2_mid==1]
timediff<-civilian$timediff[civilian$d2_mid==1]
Z<-civilian[civilian$d2_mid==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for medium offenses
m5<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m6<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2h15<-civilian$d2h15[civilian$organize==1|civilian$prepare==1|civilian$conspiracy==1|civilian$attemp==1]
timediff<-civilian$timediff[civilian$organize==1|civilian$prepare==1|civilian$conspiracy==1|civilian$attemp==1]
Z<-civilian[civilian$organize==1|civilian$prepare==1|civilian$conspiracy==1|civilian$attemp==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for medium offenses (based on crime description)
m7<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m8<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

r<-NULL
for(i in 1:8){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],3),", ",round(get(paste0("m", i))$ci[2,2],3),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],3),", ",round(get(paste0("m", i))$ci[3,2],3),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.5. Regression discontinuity of the 1956 reform on sentence length for civilian dissidents across different types of crimes------------
d2_termmm<-civilian$d2_termmm[civilian$d2_maj==1]
timediff<-civilian$timediff[civilian$d2_maj==1]
Z<-civilian[civilian$d2_maj==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for major offenses
m9<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m10<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$weapon==1|civilian$subversion==1|civilian$leakmilitary==1]
timediff<-civilian$timediff[civilian$weapon==1|civilian$subversion==1|civilian$leakmilitary==1]
Z<-civilian[civilian$weapon==1|civilian$subversion==1|civilian$leakmilitary==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for major offenses (based on crime description)
m11<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m12<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_mid==1]
timediff<-civilian$timediff[civilian$d2_mid==1]
Z<-civilian[civilian$d2_mid==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for medium offenses
m13<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m14<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$organize==1|civilian$prepare==1|civilian$conspiracy==1|civilian$attemp==1]
timediff<-civilian$timediff[civilian$organize==1|civilian$prepare==1|civilian$conspiracy==1|civilian$attemp==1]
Z<-civilian[civilian$organize==1|civilian$prepare==1|civilian$conspiracy==1|civilian$attemp==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for medium offenses (based on crime description)
m15<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m16<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_819==1]
timediff<-civilian$timediff[civilian$d2_819==1]
Z<-civilian[civilian$d2_819==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for minor offenses 
m17<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m18<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$connive==1]
timediff<-civilian$timediff[civilian$connive==1]
Z<-civilian[civilian$connive==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for minor offenses (based on crime description)
m19<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
              covs=Z,bwcheck=15,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m20<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw

r<-NULL
for(i in 9:20){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],2),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],2),", ",round(get(paste0("m", i))$ci[2,2],2),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],2),", ",round(get(paste0("m", i))$ci[3,2],2),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.6. Regression discontinuity of the 1956 reform on sentence severity for civilian dissidents across different types of offenses (alternative model specifications)------------
d2h15<-civilian$d2h15[civilian$d2_maj==1]
timediff<-civilian$timediff[civilian$d2_maj==1]
Z<-civilian[civilian$d2_maj==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for major offenses
m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m2<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m3<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m4<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m5<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m6<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2h15<-civilian$d2h15[civilian$d2_mid==1]
timediff<-civilian$timediff[civilian$d2_mid==1]
Z<-civilian[civilian$d2_mid==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for medium offenses
m7<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, bwcheck=130,
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m8<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m9<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, bwcheck=130,
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m10<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, bwcheck=130,
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m11<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m12<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2,
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

r<-NULL
for(i in 1:12){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],2),", ",round(get(paste0("m", i))$ci[2,2],2),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],2),", ",round(get(paste0("m", i))$ci[3,2],2),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.7. Regression discontinuity of the 1956 reform on sentence length for civilian dissidents across different types of crimes (alternative model specifications)------------
d2_termmm<-civilian$d2_termmm[civilian$d2_maj==1]
timediff<-civilian$timediff[civilian$d2_maj==1]
Z<-civilian[civilian$d2_maj==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for major offenses
m13<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m14<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m15<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m16<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m17<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m18<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_mid==1]
timediff<-civilian$timediff[civilian$d2_mid==1]
Z<-civilian[civilian$d2_mid==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for medium offenses
m19<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m20<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m21<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 2, bwcheck=90,
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m22<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m23<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m24<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_819==1]
timediff<-civilian$timediff[civilian$d2_819==1]
Z<-civilian[civilian$d2_819==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for minor offenses
m25<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m26<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m27<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m28<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m29<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m30<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

r<-NULL
for(i in 13:30){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],2),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],1),", ",round(get(paste0("m", i))$ci[2,2],1),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],1),", ",round(get(paste0("m", i))$ci[3,2],1),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.8. Regression discontinuity of the 1956 reform on sentence severity for civilian dissidents (donut-hole test)------------
d2h15<-civilian$d2h15[abs(civilian$timediff) >= 31]
timediff<-civilian$timediff[abs(civilian$timediff) >= 31]
Z<-civilian[abs(civilian$timediff) >= 31,
            c("age","male","islander","d2_per","weapon","subversion",
              "leakmilitary","presreyearr_l")]
m1<-rdrobust(d2h15,timediff,kernel = "tri",p = 1,covs=Z,bwselect="cerrd", bwcheck=250,masspoints = "adjust")
# 1 month donut-hole

d2h15<-civilian$d2h15[abs(civilian$timediff) >= 91]
timediff<-civilian$timediff[abs(civilian$timediff) >= 91]
Z<-civilian[abs(civilian$timediff) >=91,
            c("age","male","islander","d2_per","weapon","subversion",
              "leakmilitary","presreyearr_l")]
m2<-rdrobust(d2h15,timediff, kernel = "tri",p = 1, covs=Z,bwselect="cerrd", bwcheck=250,masspoints = "adjust")
# 3 months donut-hole

r<-NULL
for(i in 1:2){
  r<-rbind(r,cbind(round(get(paste0("m",i))$coef[1],3),
                   paste0("[",round(get(paste0("m",i))$ci[3,1],3),",",round(get(paste0("m",i))$ci[3,2],3),"]"), 
                   round(get(paste0("m",i))$bws[1,1]),
                   sum(get(paste0("m",i))$N_h)))
}
r<-data.frame(r)
r


### Table A.9. Regression discontinuity of the 1956 reform on sentence severity for civilian dissidents (with lagged dependent variable)------------
d2h15<-civilian$d2h15
timediff<-civilian$timediff
Z<-civilian[,c("age","male","islander","d2_per","weapon","subversion",
               "leakmilitary","presreyearr_l","d2h15_lm")]
#include monthly average of the dependent variable in t − 1 as a covariate

m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m2<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m3<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m4<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust",bwcheck=250) #quadratic fit, uniform kernel, CE bw
m5<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw
m6<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m7<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m8<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

r<-NULL
for(i in 1:8){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],3),", ",round(get(paste0("m", i))$ci[2,2],3),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],3),", ",round(get(paste0("m", i))$ci[3,2],3),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.10. Regression discontinuity of the 1956 reform on sentence severity for civilian dissidents (with judge fixed effects)------------
civilian <- dummy_cols(civilian, select_columns = 'd2_j1')
somejudges <- civilian %>%
  group_by(d2_j1) %>%
  filter(all(1%in% after1956)) %>%
  ungroup()
#generate dummies of judges

d2h15<-somejudges$d2h15
timediff<-somejudges$timediff
Z<-somejudges[,c("age","male","islander","d2_per","weapon","subversion",
                 "leakmilitary","presreyearr_l",
                 grep('d2_j1_', names(somejudges), value = TRUE))]
Z$d2_j1_<-NULL
Z$d2_j1_NA<-NULL
#include judge fixed effects as covariates

m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m2<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m3<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m4<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust",bwcheck=235) #quadratic fit, uniform kernel, CE bw
m5<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw
m6<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m7<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m8<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

r<-NULL
for(i in 1:8){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[2,1],3),", ",round(get(paste0("m", i))$ci[2,2],3),"]"),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],3),", ",round(get(paste0("m", i))$ci[3,2],3),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.11. Regression discontinuity of the 1956 reform on sentence severity for military cases involving generals------------
d2h15<-trial$d2h15[trial$category_e==3]
timediff<-trial$timediff[trial$category_e==3]
Z<-trial[trial$category_e==3,c("age","male","islander","d2_per","weapon","subversion",
                               "leakmilitary","presreyearr_l")]
#DV, running variable, and covariates of the RD models for generals

m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, CE bw
m2<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m3<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m4<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m5<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, triangular kernel, MSE bw
m6<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m7<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m8<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw
r<-NULL
for(i in 1:8){
  r<-cbind(r,rbind(round(get(paste0("m", i))$coef[3],3),
                   paste0("[",round(get(paste0("m", i))$ci[3,1],3),", ",round(get(paste0("m", i))$ci[3,2],3),"]"),
                   round(get(paste0("m", i))$bws[1,1]),
                   sum(get(paste0("m", i))$N_h)))
}
r<-data.frame(r)
r


### Table A.12. Difference-in-differences analyses on the effect of the 1956 reform on sentence severity------------
trial$no_military=NA
trial$no_military[trial$category_e==3&!is.na(trial$category_e)]<-0 
trial$no_military[trial$category_e==0&!is.na(trial$category_e)]<-1 
# generate a group dummy for civilian/military

shift_year <- function(year, month) {
  adjusted_year <- if (month %in% c(10, 11, 12)) year + 1 else year
  return(adjusted_year)
}
trial <- trial %>%
  mutate(year = as.numeric(format(as.Date(paste0(minym, "-01")), "%Y")),
         month = as.numeric(format(as.Date(paste0(minym, "-01")), "%m")),
         shift_year = mapply(shift_year, year, month))

gdat<-trial %>%
  group_by(no_military, shift_year, minym) %>%
  summarize(d2h15 = mean(d2h15,na.rm=T),
            age = mean(age,na.rm=T),
            male = mean(male,na.rm=T),
            islander = mean(islander,na.rm=T),
            lnd2per= mean(lnd2per,na.rm=T),
            weapon = mean(weapon,na.rm=T),
            subversion = mean(subversion,na.rm=T),
            leakmilitary = mean(leakmilitary,na.rm=T),
            presreyearr_l = mean(presreyearr_l,na.rm=T))
#aggregate all variables to yearly data, with each year running from the previous October to September of the current year

gdat<-gdat[gdat$shift_year>=1950&gdat$shift_year<=1972,] #include only data under Chiang's reign from 1950 to 1972 
gdat$after1956<-0
gdat$after1956[gdat$shift_year>=1957]<-1 #dummy of post reform
gdat$after1956<- factor(gdat$after1956, labels=c("Before 1956","After 1956"))
gdat$no_military<- factor(gdat$no_military, labels=c("Military","Civilian"))

m1<-lm(d2h15~after1956*no_military+age+male+islander+lnd2per+
         weapon+subversion+leakmilitary+presreyearr_l,data=gdat) #with group FE
m2<-lm(d2h15~after1956*no_military+age+male+islander+lnd2per+
         weapon+subversion+leakmilitary+presreyearr_l+factor(shift_year),data=gdat) #with group FE + year FE
stargazer(m1,m2)


### Figure A.1. Frequency of civilian dissident cases around Oct. 1st, 1956------------
freq<-civilian[!is.na(civilian$d2h15),]
freq<-freq[!is.na(freq$minyy),]

case<-freq %>%
  group_by(minyy,minym) %>%
  summarize(d2_all = n(),
            d2_maj=sum(d2_maj,na.rm=T),d2_mid=sum(d2_mid,na.rm=T),d2_min=sum(d2_819,na.rm=T))
#aggregate data to monthly frequencies, divided by types of offenses

case<-as.data.frame(case)
r<-names(table(case$minym))
case$minym<-factor(case$minym, levels=r)
long<-reshape(data=case,idvar="minym",varying=names(case[,3:6]),
              timevar="type",direction="long",sep="_")
long$type<-factor(long$type, labels=c("All","Major", "Medium","Minor"))

ggplot(long[case$minyy>=1955&case$minyy<=1958&!case$minym%in%c("1958-11","1958-12"),], 
       aes(minym,d2))+geom_col(col="black")+facet_grid(rows = vars(type),scales = "free", space = "free")+  
  theme( panel.grid.major =  element_line(colour = "grey90"), panel.background = element_rect(size = 1,colour = "lightgray",fill='white'),
         axis.text.x = element_text(angle = 90, vjust = 1, hjust=1.5))+
  scale_x_discrete("")+scale_y_continuous("Count",breaks = breaks_extended(8))+
  geom_vline(xintercept = 20.5, color = "#E66100")


### Figure A.2. The 1956 reform and the severity of sentences across different types of offenses (alternative model specifications)------------
d2h15<-civilian$d2h15[civilian$d2_maj==1]
timediff<-civilian$timediff[civilian$d2_maj==1]
Z<-civilian[civilian$d2_maj==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for major offenses
m1<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m2<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m3<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m4<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m5<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m6<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2, 
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2h15<-civilian$d2h15[civilian$d2_mid==1]
timediff<-civilian$timediff[civilian$d2_mid==1]
Z<-civilian[civilian$d2_mid==1,c("age","male","islander","presreyearr_l")]
#binary DV, running variable, and covariates of the RD models for medium offenses
m7<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="cerrd", p = 2, bwcheck=130,
             covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m8<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
             covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m9<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="cerrd", p = 2, bwcheck=130,
             covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m10<-rdrobust(d2h15, timediff, kernel = "tri",bwselect="mserd", p = 2, bwcheck=130,
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m11<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m12<-rdrobust(d2h15, timediff, kernel = "uni",bwselect="mserd", p = 2,
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_maj==1]
timediff<-civilian$timediff[civilian$d2_maj==1]
Z<-civilian[civilian$d2_maj==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for major offenses
m13<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m14<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m15<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m16<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m17<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m18<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_mid==1]
timediff<-civilian$timediff[civilian$d2_mid==1]
Z<-civilian[civilian$d2_mid==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for medium offenses
m19<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m20<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m21<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 2, bwcheck=90,
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m22<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m23<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m24<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

d2_termmm<-civilian$d2_termmm[civilian$d2_819==1]
timediff<-civilian$timediff[civilian$d2_819==1]
Z<-civilian[civilian$d2_819==1,c("age","male","islander","presreyearr_l")]
#continuous DV, running variable, and covariates of the RD models for minor offenses
m25<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, CE bw
m26<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, CE bw
m27<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="cerrd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, CE bw
m28<-rdrobust(d2_termmm, timediff, kernel = "tri",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, triangular kernel, MSE bw
m29<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 1, 
              covs=Z,masspoints = "adjust") #linear fit, uniform kernel, MSE bw
m30<-rdrobust(d2_termmm, timediff, kernel = "uni",bwselect="mserd", p = 2, 
              covs=Z,masspoints = "adjust") #quadratic fit, uniform kernel, MSE bw

r<-NULL
for(i in 1:30){
  r<-rbind(r,cbind(get(paste0("m", i))$coef[3],get(paste0("m", i))$ci[2,1],get(paste0("m", i))$ci[2,2],
                   get(paste0("m", i))$ci[3,1],get(paste0("m", i))$ci[3,2]))
}
r<-data.frame(r)
colnames(r)<-c("resp","bclower","bcupper","rlower","rupper")
r$group<-c(rep("DV: more than 15yr",12),rep("DV: sentence in months",18))
r$bw<-rep(c(rep("CE",3),rep("MSE",3)),5)
r$trt <-c(rep(3,3),rep(2.9,3),rep(2,3),rep(1.9,3),
          rep(3,3),rep(2.9,3),rep(2,3),rep(1.9,3),rep(1,3),rep(0.9,3))
r$ms <-rep(c("Quadratic fit, triangular kernel",
             "Linear fit, uniform kernel",
             "Quadratic fit, uniform kernel"),10)
r$group<-factor(r$group, levels=c('DV: more than 15yr',"DV: sentence in months"))

ggplot(r, aes(resp,trt,color=bw))+
  geom_pointrange(aes(xmin = bclower, xmax = bcupper,shape=bw))+
  geom_linerange(aes(y=trt-0.05,xmin = rlower, xmax = rupper),linetype=3)+
  facet_grid(vars(ms),vars(group),scales = "free")+
  theme( panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.position="top",legend.direction = "horizontal",
         axis.text.y = element_text(size =12),
         strip.text = element_text(size =12),
         panel.background = element_rect(size = 1,colour = "lightgray",fill='white'))+
  scale_y_continuous("",
                     breaks=c(2.95,1.95,0.95),
                     labels=c("Major offenses",
                              "Medium offenses",
                              "Minor offenses"))+
  scale_color_brewer("Bandwidth",palette = "Dark2")+
  scale_shape_discrete("Bandwidth")+
  geom_vline(xintercept = 0, color = "darkgray")+
  scale_x_continuous("Local average treatment effect")
#Plot the LATE point estimates and 95% CI


### Figure A.3. RD estimates with alternative bandwidths------------
d2h15<-civilian$d2h15
timediff<-civilian$timediff
Z<-civilian[,c("age","male","islander","lnd2per","weapon","subversion",
               "leakmilitary","presreyearr_l")]
#DV, running variable, and covariates of the RD models

different_h<-data.frame(h=seq(200,2500,50),
                        p=NA,l95=NA,h95=NA)
for(i in seq(200,2500,50)){
  m<-rdrobust(d2h15, timediff, kernel = "tri",h=i, p = 1, 
              covs=Z,masspoints = "adjust")
  different_h$p[different_h$h==i]<-m$coef[2]
  different_h$l95[different_h$h==i]<-m$ci[2,1]
  different_h$h95[different_h$h==i]<-m$ci[2,2]
}
#alternative bandwidths ranging from 200 to 2500

ggplot(different_h, aes(x=h, y=p)) + 
  geom_pointrange(aes(ymin=l95, ymax=h95))+
  geom_hline(yintercept=0,color="grey40")+
  theme( panel.grid.major =  element_line(colour = "grey90"), 
         panel.background = element_rect(size = 1,colour = "lightgray",fill='white'),
         legend.position="bottom")+
  scale_x_continuous("Bandwidth")+
  scale_y_continuous("RD estimate")


### Figure A.4. RD estimates for false cutoffs------------
d2h15<-civilian$d2h15
timediff<-civilian$timediff
Z<-civilian[,c("age","male","islander","lnd2per","weapon","subversion",
               "leakmilitary","presreyearr_l")]
placebo_c<-data.frame(h= seq(-2190,3285,365),
                      p=NA,l95=NA,h95=NA)
#DV, running variable, and covariates of the RD models

for(i in placebo_c$h){
  m<-rdrobust(d2h15, timediff, kernel = "tri",c=i, bwselect="cerrd",
              covs=Z,masspoints = "adjust")
  placebo_c$p[placebo_c$h==i]<-m$coef[2]
  placebo_c$l95[placebo_c$h==i]<-m$ci[2,1]
  placebo_c$h95[placebo_c$h==i]<-m$ci[2,2]
}
placebo_c$cutoff<-paste0("Oct. 1, ",seq(1950,1965,1))
#placebo tests using fake cutoffs between 1950 and 1963

ggplot(placebo_c, aes(x=cutoff, y=p)) + 
  geom_pointrange(aes(ymin=l95, ymax=h95))+
  geom_hline(yintercept=0,color="grey40")+
  theme( panel.grid.major =  element_line(colour = "grey90"), 
         panel.background = element_rect(size = 1,colour = "lightgray",fill='white'),
         legend.position="bottom",axis.text.x = element_text(angle = 90, vjust = 1, hjust=1.5))+
  scale_x_discrete("Cutoff")+
  scale_y_continuous("RD estimate")


### Figure A.5. Distribution of defendants by ideological affiliation------------
trial$communist<-0
trial$communist[grepl("匪黨",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("匪黨",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("匪黨",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("匪黨",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("匪黨",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$communist[grepl("共匪",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("共匪",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("共匪",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("共匪",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("共匪",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$communist[grepl("匪區",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("匪區",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("匪區",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("匪區",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("匪區",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$communist[grepl("匪軍",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("匪軍",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("匪軍",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("匪軍",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("匪軍",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$communist[grepl("延安",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("延安",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("延安",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("延安",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("延安",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$communist[grepl("共產黨",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("共產黨",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("共產黨",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("共產黨",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("共產黨",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$communist[grepl("毛匪",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$communist[grepl("毛匪",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$communist[grepl("毛匪",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$communist[grepl("毛匪",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$communist[grepl("毛匪",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
#Classify whether the case involves communist ideology based on court descriptions

trial$independence<-0
trial$independence[grepl("獨立",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$independence[grepl("獨立",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$independence[grepl("獨立",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$independence[grepl("獨立",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$independence[grepl("獨立",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
trial$independence[grepl("組黨",trial$d1_crime1)&!is.na(trial$d1_crime1)]<-1
trial$independence[grepl("組黨",trial$d1_crime2)&!is.na(trial$d1_crime2)]<-1
trial$independence[grepl("組黨",trial$d1_crime3)&!is.na(trial$d1_crime3)]<-1
trial$independence[grepl("組黨",trial$d1_crime4)&!is.na(trial$d1_crime4)]<-1
trial$independence[grepl("組黨",trial$d1_crime5)&!is.na(trial$d1_crime5)]<-1
#Classify whether the case involves Taiwan independence movement based on court descriptions

agg_trial<-as.data.frame(trial[trial$category_e==0,] %>%
                           group_by(minyy) %>% 
                           summarise(communist=sum(communist,na.rm = TRUE),
                                     independence=sum(independence,na.rm = TRUE)))
#aggregate data to yearly frequencies, divided by ideological affiliations

aggl <- reshape(agg_trial, 
                varying = c("communist", "independence"), 
                v.names = "count",
                timevar = "type", 
                times = c("communist", "independence"),
                direction = "long")

ggplot(aggl[aggl$minyy>=1950&aggl$minyy<=1972&!is.na(aggl$minyy),], 
       aes(x=minyy, y=count, group=type)) +
  geom_line(aes(color=type))+
  geom_point(aes(color=type))+
  scale_color_brewer(palette="Dark2",name="Dissident ideology",
                     labels=c("Communist","Taiwan independence"))+  
  theme( panel.grid.minor = element_blank(),
         panel.background = element_rect(size = 1,colour = "lightgray",fill='white'))+
  scale_x_continuous("Year")+ scale_y_continuous("Frequency")+
  geom_vline(xintercept = 1956, color = "black")
#Plot the frequencies


### Figure A.6. The estimated percentage of cases exceeding 15yr sentences------------
trial$no_military=NA
trial$no_military[trial$category_e==3&!is.na(trial$category_e)]<-0 
trial$no_military[trial$category_e==0&!is.na(trial$category_e)]<-1 
# generate a group dummy for civilian/military

shift_year <- function(year, month) {
  adjusted_year <- if (month %in% c(10, 11, 12)) year + 1 else year
  return(adjusted_year)
}
trial <- trial %>%
  mutate(year = as.numeric(format(as.Date(paste0(minym, "-01")), "%Y")),
         month = as.numeric(format(as.Date(paste0(minym, "-01")), "%m")),
         shift_year = mapply(shift_year, year, month))

gdat<-trial %>%
  group_by(no_military, shift_year, minym) %>%
  summarize(d2h15 = mean(d2h15,na.rm=T),
            age = mean(age,na.rm=T),
            male = mean(male,na.rm=T),
            islander = mean(islander,na.rm=T),
            lnd2per= mean(lnd2per,na.rm=T),
            weapon = mean(weapon,na.rm=T),
            subversion = mean(subversion,na.rm=T),
            leakmilitary = mean(leakmilitary,na.rm=T),
            presreyearr_l = mean(presreyearr_l,na.rm=T))
#aggregate all variables to yearly data, with each year running from the previous October to September of the current year

gdat<-gdat[gdat$shift_year>=1950&gdat$shift_year<=1972,] #include only data under Chiang's reign from 1950 to 1972 
gdat$after1956<-0
gdat$after1956[gdat$shift_year>=1957]<-1 #dummy of post reform
gdat$after1956<- factor(gdat$after1956, labels=c("Before 1956","After 1956"))
gdat$no_military<- factor(gdat$no_military, labels=c("Military","Civilian"))

m1<-lm(d2h15~after1956*no_military+age+male+islander+lnd2per+
         weapon+subversion+leakmilitary+presreyearr_l,data=gdat) #with group FE
m2<-lm(d2h15~after1956*no_military+age+male+islander+lnd2per+
         weapon+subversion+leakmilitary+presreyearr_l+factor(shift_year),data=gdat) #with group FE + year FE

set_theme(
  base = theme_sjplot(),
  axis.title.size = 1,
  axis.textsize.y = 1,
  legend.size = 1,
)
plot_model(m2, type = "pred", terms = c("after1956","no_military"),title="",legend.title="",
           ci.lvl =.95,colors="Dark2",
           axis.title = c("","Predicted percentage of more than 15yr sentence"))+
  ggplot2::scale_y_continuous(labels = function(x) insight::format_value(x,
                                                                         as_percent=TRUE,
                                                                         digits=0))
#Plot model estimates


### Figure A.7. The dynamic treatment effect of the 1956 reform on sentence severity------------
gdat$no_military<- as.numeric(gdat$no_military)
est_did <-feols(d2h15 ~age+male+islander+lnd2per+
                  weapon+subversion+leakmilitary+presreyearr_l+
                  i(shift_year, no_military, 1956) |no_military+shift_year, 
                gdat) #dynamic model with treatment and year fixed effects
iplot(est_did, ci.lwd = 2, pt.cex = 2, main = "", 
      xlab = "Year",  ref.line.par = list(lwd = 2))


### Figure A.8. Relative frequency of presidential reviews for civilian cases------------
overrule<-civilian[!is.na(civilian$d2h15),]
overrule$presn0<-overrule$presn
overrule$presn0[overrule$pres==0]<-0
overrule$presn0[overrule$presn%in%c(2,3,4,5,6)]<-2
#Indicate whether the President reviewed the case once or multiple times

agg_trial<-as.data.frame(overrule %>%
                           group_by(minyy,presn0) %>% 
                           summarise(ncase=n())%>%
                           mutate(rel_freq = ncase / sum(ncase),
                                  freq = ncase))
agg_trial$presn0<-factor(agg_trial$presn0)
#Aggregate data to yearly frequencies, divided by the number of presidential reviews

ggplot(agg_trial[agg_trial$minyy>=1950&agg_trial$minyy<=1972&agg_trial$presn0%in%c(1,2),], 
       aes(minyy,rel_freq))+
  geom_col(aes(fill = presn0),colour="white",linewidth=0.1)+
  scale_fill_brewer(palette = "Set2",
                    name="Presidential review",labels=c("Review once","Multiple reviews"))+  
  theme(panel.grid.major = element_line(size = 0.5,colour = "gray90"),
        panel.background = element_rect(size = 1,colour = "lightgray",fill='white'))+
  scale_x_continuous("Year")+ scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                                                 "Proportion")+
  geom_vline(xintercept = 1956.5, color = "black")
#Plot the frequencies


### Figure A.9. Frequency of civilian dissident cases (1950–1972)------------
agg_trialall<-as.data.frame(civilian %>%
                              group_by(minyy) %>% 
                              summarise(ncase=n()))
#Aggregate civilian data to yearly frequencies

ggplot(agg_trialall[agg_trialall$minyy>=1950&agg_trialall$minyy<=1972,], 
       aes(minyy,ncase))+geom_col()+
  theme(panel.grid.major = element_line(colour = "lightgray",
                                        size = 0.5),
        panel.grid.minor = element_line(colour = "lightgray",
                                        size = 0.2),
        panel.background = element_rect(size = 1,colour = "lightgray",fill='white'))+
  scale_x_continuous("Year")+ scale_y_continuous("Frequency")+
  geom_vline(xintercept = 1956.5, color = "#E66100")
#Plot the frequencies
